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ABSTRACT 

Boussinesq-type equations for weakly nonlinear, weakly dispersive waves 
have been used extensively to model wave shoaling on beaches. Deterministic 
Boussinesq models cast in the form of coupled evolution equations for the 
amplitudes and phases of discrete Fourier modes (Freilich and Guza, 1984) describe 
the shoaling process accurately for arbitrary incident wave conditions, but are 
numerically cumbersome for predicting the shoaling evolution of continuous spectra 
of natural wind-generated waves. Here an alternative stochastic formulation of a 
Boussinesq model (Herbers and Burton, 1996, based on the closure hypothesis that 
phase coupling between quartets of wave components is weak) is implemented that 
predicts the evolution of a continuous frequency spectmm and bispectrum of waves 
normally incident on a gently sloping beach with straight and parallel depth 
contours. The general characteristics of the model are examined with numerical 
simulations for a wide range of incident wave conditions and bottom profiles. 
Stochastic and deterministic Boussinesq model predictions are compared to field 
observations from a cross-shore transect of bottom pressure sensors deployed on 
a barred beach near Duck, NC, during the recent DUCK94 Experiment. Predictions 
of the two models are similar and describe accurately the observed nonlinear 
shoaling transformation of wave spectra. 
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I. INTRODUCTION 



Wind-generated surface gravity waves are the principal driving force of nearshore 
fluid motions (e.g., longshore currents, rip currents, and undertow) and sediment transport 
(e.g., erosion and accretion of beaches, and the formation of bars and cusps). As waves 
shoal onto beaches, amphmdes increase, wavelengths decrease, and propagation directions 
refract towards normal incidence. These linear propagation effects are readily observed 
and well understood. Additionally, pronounced nonlinear effects in shallow water cause 
a dramatic transformation of wave shapes from initially smooth, nearly sinusoidal profiles, 
to asymmetric, pitched forward profiles characteristic of near-breaking waves. The 
mechanism for this transformation is nonlinear triad interactions in which two primary 
wave components with frequencies and (O2 excite a secondary wave component with 
the sum (co 2+0)2) or difference ((02-0)2) frequency. The nonlinear energy transfers in 
these interactions not only broaden the spectrum but also change the statistical properties 
of the waves. Whereas the primary wave components incident from deep water are 
approximately statistically independent, the newly formed secondary components are 
"phase locked" to the primary waves that excite them. Even relatively weak secondary 
components significantly change the shapes of waves in shallow water. While the incident 
waves and the nonlinearly excited higher frequency waves are predominantly dissipated 
in the surf zone, the nonlinearly excited lower frequency (infragravity) wave components 
reflect from the beach and often dominate wave runup at the shoreline. Accurate 
prediction of the nonlinear shoaling transformation of waves is critical to both naval (e.g.. 
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amphibious landings and mine warfare) and civilian (e.g., mitigation of beach erosion) 
operations in the littoral zone. 

Li deep water (Kh » 1, where k is the wavenumber and h is the water depth) and 
intermediate depths {Kh = 0(1)) triad interactions are non-resonant. The nonlinearly 
excited secondary waves remain small ("bound" waves) and are well described by fini te 
depth theory based on Stokes perturbation expansion for small wave steepness (Phillips 
1960; Hasselmann 1962; Herbers et al., 1992,1994; and many others). In shallow water 
{kh < 1), where triad interactions are near-resonant, finite depth theory is valid only for 
small values of the Ursell number, U,. = a/i^h^ (Ursell, 1953) where a is the wave 
amplitude, and this condition is usually violated on natural beaches. 

Models for waves in shallow water are usually based on the Boussinesq equations 
which assume that a/h (nonlinearity) and (x/t)^ (dispersion) are both small and of the 
same order (i.e., =0(1)). These equations are surprisingly robust, even for large t/^ 

values typically observed in near-breaking waves. Peregrine (1967) extended the 
Boussinesq equations to varying depth, and these equations form the basis of most wave 
shoaling models. Freilich and Guza (1984) developed a frequency domain (i.e., neglecting 
directional spreading effects) Boussinesq model that predicts accurately the energy 
transfers to higher ft'equencies and associated wave shape changes on natural beaches 
(Elgar and Guza, 1985a; Elgar et al., 1990a). This deterministic model is initialized at an 
offshore boundary by a discrete Fourier representation of the incident waves. A coupled 
set of evolution equations for the amplitudes and phases of the Fourier modes is solved 
numerically to evaluate the shoaling transformation of the wave train. The results of these 
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integrations are subsequently averaged over many realizations with random initial 
amplitudes and phases to obtain spectral statistics. Although the shoaling evolution of 
wave spectra is accurately reproduced, this method is cumbersome for practical 
applications, requiring large computing resources and a detailed specification of incident 
wave conditions at the offshore botmdary that is often not available. In addition, the 
extension of this approach to include directional spreading effects (i.e., two dimensions) 
is far from straightforward owing to the large number of modes required in discrete mode 
simulations of continuous frequency-directional wave spectra. 

The initialization of wave shoaling models requires local offshore measurements 
(e.g., wave followmg buoys) or predictions from a global or regional wind-wave 
generation model. These larger scale models (e.g., WAM; The WAMDI group, 1988) are 
typically stochastic formulations based on an energy balance, and predict wave spectra 
rather than individual wave profiles. Routine measurement systems provide similar 
statistical information. 

Recently, Berbers and Burton (1996) derived an alternative stochastic formulation 
of Boussmesq shoaling evolution equations, for directionally spread waves propagating 
over a gently sloping beach with straight and parallel contours. Based on the closure 
hypothesis that phase-coupling between quartets of wave components is weak, Berbers 
and Bmton derived a coupled set of evolution equations for the wave spectrum and 
bispectrum. The bispectrum describes the degree of coupling and the phase-relationship 
in triads of nonlinearly interacting wave components (Basselmann et. al., 1963), and is 
used extensively in other nonlinear process studies (e.g., economics, brain-wave emission. 
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and plasma physics). In deep water and intermediate depths, the bispectrum is completely 
determined by the local spectrum, and enables the detection of relatively weak phase- 
coupled, forced secondary waves that are concealed in the spectrum by more energetic 
freely propagating primary waves (e.g., Herbers et al., 1992, 1994). In shallow water, the 
bispectrum evolves strongly and describes, in a statistical sense, the shape evolution of 
shoaling waves (e.g., Elgar and Guza, 1985b; Elgar et al., 1990a). 

Here, a numerical implementation of a one-dimensional stochastic Boussinesq 
model is presented in which directional spreading effects are neglected. This formulation 
allows for simple illustration of stochastic model characteristics and comparisons to field 
data and existing one-dimensional models. Energy transfers to higher frequencies in sum 
triad interactions are insensitive to directional spreading angles of incident waves (e.g., 
Herbers and Burton, 1996), and thus can be predicted accurately with a one-dimensional 
model. However, energy transfers to infragravity frequencies in difference triad 
interactions are significantly reduced for large directional spreading angles (Herbers and 
Burton, 1996). Additionally, the reflection from shore and refractive trapping of 
infragravity waves (Herbers et al., 1995) is neglected. Hence, infragravity waves are only 
crudely represented in the present model formulation. 

Following a review of the stochastic formulation of Boussinesq wave shoaling 
equations, the numerical model implementation is described in Chapter II. The 
dependence of wave shoaling evolution on the nonlinearity, spectral shape of incident 
waves and the beach profile is examined through numerical simulations in Chapter III. 
Stochastic and deterministic (the Freilich and Guza model) Boussinesq predictions are 
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compared to field data collected on a natural beach near Duck, NC in Chapter IV, 
followed by conclusions in Chapter V. 
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II. A STOCHASTIC BOUSSINESQ MODEL 



HerbeK and Burton (1996) derived a stochastic formulation of Boussinesq wave 
shoaling equations for directionally spread waves propagating over a beach with straight 
and parallel depth contours. Under the third-order closure hypothesis that phase-coupling 
between quartets of wave components is weak, the statistical properties of the waves are 
described by a coupled set of evolution equations for the frequency (6>)-alongshore 
wavenumber (/) spectrum E(o),l) and bispectrum If directional 

spreading is neglected (i.e., I = 0), these equations (22a and 22b in Herbers and Burton, 
1996) reduce to: 






^ / / 

— 5(o) ,G)-a) ) 
dx 



3 dh 
Ah dx 






( 2 ) 



- i- 






(o'e{(o -( o')E((n)+icd -u>')E(<x)')Ei(ji)-(A)Eioi^)Ei(o -oi') 



where the x-axis points onshore, E((o) and are the frequency spectnun and 

bispectrum, h(x) is the water depth, g is gravity, and IM{ } indicates the imaginary part. 
The integrals of the (density) spectrum E and bispectrum B over all frequencies yield the 
mean square <t| > and mean cube <ti > of the surface elevation ri(x,r): 
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<T|^> = J E((ji)d(A> 



<n > - 



The first term on the right-hand side of equations (1) and (2) represents linear 
shoaling effects. The nonlinear transfers in the energy spectrum are controlled by the 
imaginary part of the bispectrum (the integral on the right-hand side of equation (1)). The 
energy product terms in equation (2) represent the changes in the imaginary part of the 
bispectrum due to the three possible nonlinear interactions (one sum interaction and two 
difference interactions) within the (o)',co-co',co) triad. The second term on the right-hand 
side of equation (2) represents the detuning of the interactions from resonance (i.e., 
changes in the phase of the bispectrum owing to weak dispersion). In the limit of small 
amplitudes and bottom slope, steady solutions for E(co) and smoothly match 

the second-order bound-wave solutions of dispersive finite depth theory (Herbers and 
Burton, 1996). 

The model is initialized with a spectrum E{o)) specified at the offshore boimdary 
of the model domain (e.g., from nearby measurements or a regional model prediction) and 
the corresponding second-order finite depth theory expression for the bispectrum 
0 )'). 
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( 3 ) 



= 2[D((i)^,(i)-Ci)')£(a)^)jE(<o-(i)^) 

+ D((o',-(o)E((i>')E((i)) + D((o-(o^,-(i))E(o)-o)')E((o)] 

The coupling coefficient D is given by 

(g) + (x> f 

D(G)j,a)p = 5 

g(K^ + K2)tanh[(Kj + K^)h] - ((o^ + 0)^ 
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where his the local water depth, and the frequencies (Oj, o >2 aud wavenumbers Kj, K 2 (?q 
< 0 for £ 0 , < 0) of the interacting primary wave components obey the linear theory 
dispersion relation co, = gK-tanh(K;^)(Hasselmann, 1962; Hasselmann et al.,1963). 
Using the symmetry relations (Hasselmann et al., 1963): 

£((o) = £(-(o) 

5((0^,(i)-(0^) = 5((i>- (0^,(0^) = 5 *(-(0^,(0 ^-(o) 

= 5(g) ^,-(0) = 5(g)-g)^,-(i)) 



where * indicates the complex conjugate, the integral term in equation (1) can be 
expressed as the sum of two integrals over positive frequencies 



9 



J‘lM{B(o)',(o-(o)}d(o' = J'lM{B(o)[o)-(o)}d(o' - 2j'lM{B(o)',(o)}d(o' (5) 



that represent the energy transfers to frequency 6) resulting from the sum interaction of 
frequencies co', co-o)', and the difference interaction of frequencies 0)+co', co', 
respectively. Hence integrations of the spectrum and bispectrum evolution equations (1), 
(2) can be restricted to positive frequencies (o)', o)-co', co > 0). 

The spectrum and bispectrum are discretized: 

(0 = nAo) for n = 1, 2, N 

n 

E = £(o) ) for n = 1, 2, N 

n n 

B = R + il = B((j3 ,(o ) for n = 1, 2, N-1 and m = 1, 2, ... N-i 

nm nm nm n m 

where A O) is the bandwidth and 6)^ the highest frequency included in the computations. 
With these definitions, equations (1) and (2) reduce to a linear set of ordinary 
differential equations which can be written in the general form 

dY 

— = F(Y) 
dx 

where the elements of y the discretized spectrum (£„) and bispectrum 
values and F(T) incorporates the corresponding right-hand side of equations (1) and (2). 
This system of equations is solved using the Bulirsch-Stoer method, a variant of 
Richardson extrapolation to the limit that uses adaptive stepsize control (Press et al., 
1992; Bulirsch and Stoer, 1966). 
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III. SIMULATIONS 



To investigate the general model characteristics and the dependence of wave 
shoaling evolution on the nonlinearity, spectral shape and bottom profile, extensive 
numerical simulations were carried out with simple incident wave spectra of the general 
form: 



E^(f) 



j,n 

E—txp 

4 



( 

n I / 



1-n 









'lY" 

v4j 






n 

E sech 

K 



n 



f 



( 6 ) 



(7) 



where Ej and Ej are single-sided spectra {E(j) = 47t£((o) with /= o)/27t), E is the surface 
elevation variance <t| > and the parameter n defines the width of the spectrum. Ej is a 
broad spectrum with a power law high-frequency tail characteristic of actively generated 
wind waves. £2 has a narrow symmetric exponential shape, characteristic of remotely 
generated swell. All model simulations were initialized in a depth h = 6m with a spectral 
peak frequency^ = 0.07 Hz. The corresponding wave number Kp = 0.058 m'^ yields a 
representative value of the dispersion parameter Kh = 0.35 at the offshore boxmdary. 
Example simulations of the shoaling of a broad sea spectrum (Ej with n = 5, the Pierson- 
Moskowitz [1964] spectral shape) and a narrow swell spectrum (£2 with n = 20, the full 
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width at half maximum power is 0.009 Hz) are shown in Figures 1-5 for different bottom 
profiles and incident significant wave heights The initial bispectrum was 

obtained by substituting the initial spectrum in the finite depth theory relation (equation 

3). 

The shoaling evolution of a narrow swell spectrum with significant wave heights 
of 0.05 and 0.5 m on plane beaches with slopes of 1:30 and 1:300 is shown in Figure 1. 
All four simulations show the familiar growth of harmonic peaks at frequencies 2fp, 3fp, 
••• and an infragravity peak at about 0.01 Hz. Even for the small = 0.05 m waves, 
harmonic spectral levels are significant (up to 10 % of the primary peak level) in 1.5 m 
depth. Although the nonlinearity remains weak ([2£]^^/ h, a representative value of a/h, 
varies between 0.003 and 0.017), the Ursell number (U^ = [2£]^^^/ Kp^h^) increases from 
0.024 in 6 m depth to a relatively large value of 0.58 in 1.5 m depth. As expected, the 
shoaling evolution is much stronger for the larger = 0.5 m waves with harmonic 
spectral levels that are comparable to the primary peak levels in 1.5 m depth. In these 
simulations alh increases from 0.03 in 6 m depth to 0.17 in 1.5 m depth, and Uj. increases 
from 0.24 to 5.8. In both the = 0.05 and 0.5 m cases, stronger growth of harmonic and 
infragravity peaks is predicted on a gentle 1:300 slope than on a steep 1:30 slope. 
Eventually (Figure 1 d, f) nonlinear energy transfers fill the valleys between harmonic 
peaks and the spectrum flattens, similar to simulations with a deterministic Boussinesq 
model reported by Elgar et al. (1990b). 

The energy transfers and their dependence on the bottom slope are further 
illustrated in Figure 2 with normalized bispectrum predictions in 2 m depth. The 
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normalized bispectrum: 



w,./p = 









( 8 ) 



o 

with = 8 ti 6 ) 2 ), defined analogous to E(j), indicates a relative measure of 

phase coupling between frequencies /;, and /; + ^2 (Berbers et al., 1992). Positive 
values of the imaginary part of the bispectrum indicate energy transfers to higher 
frequencies through sum interactions, whereas negative values indicate energy transfers 
to lower frequencies through difference interactions (equations 1, 5). All four simulations 
show strong coupling at (fj, = (0.07, 0.07) Hz associated with the fp, fp, 2fp sum 
interaction and at (0.07, 0.14) Hz (the 7 ^, 2fp, 3fp sum interaction). The larger wave and 
gentle bottom slope simulations also show coupling to higher harmonics (e.g., the (0.14, 
0.14), (0.21, 0.07), and (0.21, 0.14) Hz peaks). Whereas the imaginary part of b is small 
on the gentle 1:3(X) slope (i.e., peaked but nearly symmetric wave shapes, Elgar and 
Guza, 1985) the real part of b is small on the steep 1:30 slope (i.e., pitched forward wave 
shapes). The energy transfer rate (proportional to the imaginary part of the bispectrum, 
equation 1 ) is smaller on the gentle slope than on the steep slope but the longer 
interaction distances cause stronger cumulative spectral evolution (Figure 1). 

Simulation results of the shoaling of a broad spectrum with the same initial 
significant wave heights (0.05 and 0.5 m) and beach slopes (1:30 and 1:300) are shown 
in Figure 3. The spectral evolution is much weaker than in the narrow spectra simulations 
because the principal effect of triad interactions is to spread energy to frequencies where 
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spectral levels are relatively low. In the = 0.05 m simulations (Figure 3 a, c, e) the 
nearly uniform increase in spectral levels at frequencies ^ 0.05 Hz is a linear shoaling 

(decrease in wave group speed) effect ( the first term on the right-hand side of equation 
1). The 1:300 slope simulation shows slightly larger growth of spectral levels above about 
2fp that is the result of sum interactions. The larger wave = 0.5 m) simulations 
(Figure 3 b,d,f) show the expected stronger nonlinear evolution. Although the spectrum 
remains featureless, nonlinear interactions cause a flattening to a nearly white spectrum 
in 1.5 m depth, similar to the narrow spectrum simulations (Figure If). The dependence 
on bottom slope is quaUtatively similar to the narrow spectrum results with larger 
cumulative energy transfers on a gently sloping beach (compare the solid and dashed 
curves in Figure 3). 

The shoaling evolution of a broad spectrum of waves (initial = 0.5 m) over 
three different bottom profiles is compared in Figure 4. All three profiles start with a 
gently sloping (1:300) section from 6 to 3 m depth to let the waves evolve to a shallow 
water regime with significant nonlinear energy transfers. From 3 m depth shoreward, the 
waves either continue to shoal on a 1:300 slope to 1 m depth ("beach" case), propagate 
the same 600 m distance in (constant) 3 m depth ("flat" case), or unshoal over a -1:200 
section back to 6 m depth ("bar" case). In contrast to the beach case, the spectral 
evolution between x = 900 and x = 1500 m is weak over the flat bottom, and on the back 
side of the bar high frequency spectral levels are reduced to approximately the initial 
levels in 6 m depth. At infragravity frequencies spectral levels continue to increase on all 
three profiles, but the growth is strongest on the beach profile and weakest on the bar 
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profile (Figure 4). 

Bispectra after 150 m of evolution on the three different profiles are compared in 
Figure 5. The predicted small imaginary part of the bispectrum remains positive on the 
beach profile, allowing for continued energy transfers to high frequencies (equations 1, 
5). On the flat bottom where the spectral evolution is weak, the imaginary part of the 
bispectrum shows small alternating positive and negative peaks. On the negative slope 
section of the bar profile, positive values of the imaginary part of the bispectrum (i.e., 
energy transfers to higher frequencies) evolve to negative values (i.e., energy transfers to 
lower frequencies) over a wide range of frequencies, causing a reversal in nonlinear 
energy transfers towards lower frequencies (Figure 4). 

In simulations of waves propagating over a flat bottom and into deeper water 
(Figure 4), small undulations are noted in the spectra that grow with distance and 
eventually develop into instabilities. A sensitivity analysis of the numerical solutions to 
variations in frequency bandwidth, the error tolerance of the numerical integration 
routine, the maximum frequency, and using different extrapolation techniques (i.e., 
polynomial and rational extrapolation), yielded identical features in all calculations. 
Energy was also conserved in these simulations to a high degree of accuracy. These 
numerical tests indicate that the predicted growing undulations in the spectrum are true 
features of the spectral and bispectral evolution equations and not caused by numerical 
truncation errors. However, the Boussinesq equations, truncated at second-order in 
nonlinearity and thus valid only over Oialh)'^ wavelengths (Freilich and Guza, 1984), do 
not describe accurately the longterm evolution of these moderately energetic waves. 
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Hence, the undulations in the spectrum may not be physically real but possibly result 
from the breakdown of the weakly nonlinear approximation used in the present model. 
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IV. COMPARISONS TO FIELD OBSERVATIONS 



Extensive field observations of wave shoaling were obtained during the fall of 
1994 in the DUCK94 nearshore field experiment (Elgar et al., 1996). A cross-shore 
transect of 15 SPUV systems, each consisting of a co-located pressure transducer, 
bidirectional electromagnetic current meter, and sonar altimeter were deployed on a sandy, 
barred beach near Duck, North Carolina. The 350 m long transect extended from the 
shoreline to about 6 meters depth (Figure 6). The sample frequency of all instruments was 
2 Hz. Sea surface elevation spectra with approximately 120 degrees of freedom were 
estimated from three-hour-long pressure records using a linear theory depth correction. 

The present analysis of four case studies is focused on benign wave conditions 
(incident wave significant heights ranged from 0.4 - 0.8 m) when the surf zone was 
confined to the beach face at the shoreward end of the transect. These observations span 
a 2 week period in September with small bathymetric changes. Differences between the 
depth profiles of the case studies (Figure 6) are primarily due to tidal sea level 
fluctuations. The beach profile is characterized by a sandbar located about 120-140 m 
from the shoreline and submerged approximately 2.2-2.5 m below the mean sea siuface. 
The bottom slope is approximately 1:80 seaward of the sandbar. Shoreward of the 
sandbar, the profile deepens slightly (20-40 cm) into a relatively flat, 80 m wide trough 
that extends to the steep (1:10) beach face. The beach profile used in the Boussinesq 
model computations was obtained through linear interpolation of the depth estimates 
extracted from the sonar altimeters (Gallagher et al., 1997). 
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Stochastic Boussinesq model predictions of wave spectrum evolution in the four 
case studies are compared to the observed spectrum evolution and deterministic 
Boussinesq model predictions (provided by Dr. Steve Elgar) in Figures 7, 8, 12 and 15. 
The stochastic model predictions for 15, 21 and 24 September were initialized with the 
spectrum measured at the furthest offshore pressure sensor, pl9. The 10 September case 
was initialized with pressme sensor pl8 because sensor pl9 malfunctioned during this 
run. The deterministic model predictions were initialized with the measured pressure time 
series (see Elgar et al. 1996 for further details) at the same locations. 

In all four cases (and other case studies not shown) predictions of both models 
are in excellent agreement with the observed wave shoaling evolution. The narrow swell 
spectrum cases (10 and 15 September) show the familiar amplification of harmonic peaks. 
On 10 September the incident wave ^ectrum is dominated by a narrow swell peak (peak 
frequency = 0.075 Hz), with a broader, but relatively small sea peak at 0.12 Hz (Figure 
7). Energy is transferred from the swell peak frequency 7^ to higher frequencies through 
sum triad interactions resulting in distinct harmonic peaks at 2fp (0.15 Hz; driven by fp, 
fp interactions), 3fp (0.23 Hz; fp, 2fp interactions) and 4^ (0.3 Hz; fp, 3fp and 2fp, 2fp 
interactions). Close to shore the small 0.12 Hz incident sea peak is completely submerged 
in the 2fp swell harmonic (Figure If). 

On September 15 the incident wave spectrum is distinctly bimodal with a narrow 
swell peak fp ~ 0.06 Hz) and a slightly smaller sea peak at twice the swell frequency (2fp 
-0.12 Hz) (Figure 8). As these wave systems propagate over the shallow sandbar, large 
nonlinear energy transfers in sum interactions yield clearly distinguishable harmoitic peaks 
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at 3fp, 4fp and 5fp (compare Figure 8a and 8c). The cross-shore evolution of spectral levels 
at frequencies fp, 2fp, 3fp, 4fp and 5fp is shown in Figure 9. Partial reflection of the 
dominant 0.06 Hz swells from the beach is evident in the large cross-shore energy 
variations (i.e., standing wave patterns) observed close to shore. Good agreement between 
the observed and predicted growth of higher-frequency harmonics indicates that nonlinear 
energy transfers are insensitive to weak reflections from shore (Elgar et al., 1996). 

Both the observed and predicted bispectra on September 15 show the expected 
shoaling transition from real values (i.e., skewed but symmetric wave profiles) (Figure 
10a and 10b) to imaginary values (i.e., pitched-forward wave profiles) (Figure 10c and 
lOd). Although the observed and predicted bispectral levels generally agree within the 
considerable uncertainty in estimates extracted from short field data records (records 
longer than 3 hours could not be used owing to tide-induced depth changes), they differ 
in detail at the shallower sites (Figure 1 1). In frequency pairs involving the 0.06 Hz swell 
peak, the observed bispectrum shows a dramatic shift from real to imaginary values 
between sensors p3 and p4 (separated by only 15 m) that is absent in the model results. 
This biphase shift is caused by the partial reflection of the 0.06 Hz swell from shore 
(Figure 9) that is not incorporated in the model predictions . Midway between nodes and 
antinodes (e.g., sensor p4) the reflected components are 90° out of phase with the incident 
components, causing large biphase shifts in triads involving the standing wave component. 

The September 24 spectrum, is broader and more energetic, with a peak frequency, 
^ = 0.1 Hz (Figure 12). Sum interactions transfer energy to a broad range of higher 
frequencies, causing a broadening of the spectrum rather than the development of distinct 
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harmonic peaks observed on September 10 and 15. The observed and predicted spectral 
levels at 2>fp decrease sharply between the bar crest (x = 240 m) and the slightly deeper 
trough (x = 300 m), whereas energy levels at 2fp continue to increase (Figure 13). These 
results suggest that energy is transferred back to lower frequencies as waves travel over 
the sandbar into deeper water similar to the simulation results (Figures 4, 5). Observed 
and predicted bispectra (Figure 14) show a clear transition from positive imaginary parts 
seaward of the bar crest to negative imaginary parts shoreward of the bar crest at high 
frequencies that is consistent with this reversal in the nonlinear energy transfer (equations 
1, 5). 

In contrast to the September 10, 15 and 24 case studies, the shoaling evolution of 
the broad, featureless spectrum observed on September 21 (Figure 15) is weak even 
though the nonlinearity is comparatively strong in this case (the incident wave significant 
height is 0.8 m). Sum and difference interactions cause significant energy transfers, but 
do not strongly affect the spectrum because the interactions, spread over a wide frequency 
range, tend to cancel out in a broad spectrum. Predicted bispectral levels (not shown) are 
low, consistent with the observations. 

Discrepancies between the stochastic model predictions and the observations, and 
between the deterministic model predictions and the observations are roughly comparable 
(observed and predicted spectral levels generally agree within about a factor of four) but 
differ in detail. Initially the deterministic model tends to overpredict energy transfers to 
higher frequencies whereas the stochastic model predictions are in close agreement with 
the observed spectra (e.g.. Figures 7a-c, 8a-b). During the later stages of shoaling 
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evolution, the stochastic model tends to overpredict high-frequency spectral levels, 
whereas the deterministic model predictions are close to, or in some cases, underpredict 
(e.g.. Figure 12 c, d) the observed spectral levels. Some of these differences may be the 
result of the different formulation of dispersion effects in the two models. The stochastic 
model uses Peregrine’s (1967) "consistent" approximation, whereas the deterministic 
model uses Freilich and Guza’s (1984) "dispersive" approximation. Comparisons of 
deterministic Boussinesq models to field data reported by Freilich and Guza (1984) show 
similar trends of overprediction (for the "consistent" approximation) and underprediction 
("dispersive" approximation) of high frequency spectral levels. Other possible explanations 
for small differences between the deterministic and stochastic model predictions are the 
different way the models initialize third-order statistics (the stochastic model uses second- 
order finite depth theory whereas the deterministic model uses measured time series) and 
the statistical closure of the stochastic model. Additionally, dissipation (neglected in both 
models) and higher-order nonlinear effects likely contribute significant errors in the 
predictions close to shore. 

The predicted shoaling amplification of low frequency (<0.06 Hz) spectral levels 
is generally in reasonable agreement with the observations, even though the model is 
obviously inadequate at infragravity firequencies. It is well known that reflection from 
shore (e.g., Elgar et al., 1994) and refractive trapping in deeper water (e.g., Huntley et al., 
1981; Herbers et al., 1995 a, b) are important at infragravity frequencies, and these effects 
are not incorporated in the models. Furthermore, the energy transfers to infragravity 
frequencies are sensitive to directional spreading effects that are neglected here (Herbers 
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et al., 1995; Herbers and Burton, 1996). Nevertheless, the roughly comparable observed 
and predicted infragravity energy levels suggest that nonlinear triad interactions is a 
plausible mechanism for the transfers of energy to infragravity frequencies in shallow 
water. 
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V. SUMMARY AND CONCLUSIONS 



A stochastic model for the shoaling of waves on a beach with straight and parallel 
depth contours is presented, based on a third-order closure of Boussinesq equations 
(Herbers and Burton, 1996). The model includes nonlinear triad interactions in which two 
primary wave components with frequencies 0)j and (O 2 , excite a secondary wave 
component with the sum or difference (o)j- 0 ) 2 ) frequency. Neglecting directional 

spreading effects, a coupled set of evolution equations for the wave spectrum and 
bispectrum is solved with standard numerical integration techniques. The model is 
numerically efficient and requires only an estimate of the incident wave spectrum for 
initialization that is often readily available from a nearby buoy or a regional wave model 
prediction. The bispectrum is initialized with a local prediction based on second-order 
finite depth theory. 

Extensive numerical simulations were performed to examine the model 
characteristics and the dependence of wave shoaling evolution on nonlinearity, spectral 
shape and bottom profile. In simulations with strong nonlinearity, both narrow and broad 
spectra tend to evolve to a flat featureless spectnun (Figures If, 3f). Simulations of 
narrow spectra with peak frequency show the familiar growth of harmonic peaks at Ifp, 
2>fp, ••-. In simulations with broad spectra comparable energy transfers to higher 
frequencies occur, but since the interactions are spread over a wide frequency range, the 
spectra remain featureless (Figure 3). On gently sloping beaches predicted nonlinear 
energy transfer rates are weaker than on steep beaches but cause stronger cumulative 
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spectral evolution (Figures 1,3)- On steep slopes, predicted bispectra have relatively large 
imaginary parts characteristic of pitched forward wave shapes, whereas the predominantly 
real bispectral values predicted on gentle slopes indicate symmetric wave profiles . These 
characteristics are qualitatively consistent with wave shape evolution observed prior to 
breaking on natural beaches. Simulations of waves propagating over a bar into deeper 
water show decreasing high-frequency spectral levels (Figure 4). The predicted bispectra 
indicate a reversal in nonlinear energy transfers with difference triad interactions 
transferring energy back towards lower frequencies. 

Stochastic and deterministic (Freilich and Guza, 1984) Boussinesq model 
predictions were compared to extensive field observations of wave shoaling on a natural 
barred beach. Although predictions of the two models differ in detail, the overall 
agreement with the observed wave spectrum evolution is comparable. Both models predict 
accurately the nonlinear transfer of energy to higher frequencies for a range of incident 
wave conditions (Figures 7, 8, 12, 15). These results are similar to earlier studies using 
deterministic Boussinesq models on nearly plane California beaches (Freilich and Guza, 
1984; Elgar and Guza, 1985). Although spectral levels at high frequencies generally 
increase as waves propagate shoreward owing to sum triad interactions, in one case a 
decrease in high-frequency spectral levels was observed shoreward of the sandbar, 
consistent with difference interactions predicted by both models (Figures 12, 13). These 
observations support the simulation result that energy can be transferred back to incident 
wave frequencies in regions of gradually increasing depth. 

In conclusion, a numerically efficient stochastic Boussinesq model was developed 
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that predicts the evolution of wave spectra and bispectra on beaches. Comparisons of 
model predictions to extensive field observations firom a natural beach show excellent 
agreement. In the future, the stochastic wave shoaling model will be extended to 
directionally spread waves based on similar evolution equations for the frequency- 
alongshore wavenumber spectrum and bispectrum (Herbers and Burton, 1996). 
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APPENDIX 



Figure 1. Niimerical simulations of the shoaling evolution of a narrow spectrum 
of waves (equation 7 with n = 20) over a plane beach. The model is initialized in 6 m 
depth. Results are shown in 4 m (upper panels), 2 m (middle panels), and 1.5 m (lower 
panels) depth for incident wave significant heights of 0.05 m (left panels) and 0.5 m 
(right panels) and beach slopes of 1:300 (solid lines) and 1:30 (dashed lines). The initial 
spectrum is indicated in each panel with a dotted line. 

Figure 2. Normalized bispectra bifj, (equation 8, units Hz'^^^) predicted in 2 
m depth in the simulations described in Figure 1. The real and imaginary parts of b are 
shown in the lower and left quadrants, respectively. Contour levels are: ± 1, 3, 5. 

Figure 3. Numerical simulation of the shoaling evolution of a broad spectrum of 
waves (equation 6 with n = 5) over a plane beach (same format as Figure 1). 

Figure 4. Simulation of the evolution of a broad spectrum (equation 6 with n = 
5), with an initial = 0.5 m in 6 m depth, over 3 different bottom profiles: a beach (left 
panels), a flat section (center panels), and a bar (right panels). Predictions are shown at 
X = 900, 1050, and 1200 m (indicated by asterisks in the upper panels). The dotted line 
indicates the initial spectrum at ;c = 0. 

Figure 5. Normalized bispectra predicted in the simulations described in Figure 
4. (a) X = 900 m (same for all 3 profiles), (b) x = 1050 m on the beach profile, (c) x = 
1050 m on the flat section, (d) x = 1050 m on the bar profile. The format of the panels 
is the same as in Figure 2. 
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Figure 6. Depth profiles and sensor locations of the 4 field data case studies. 

Figure 7. Comparison of observed (solid line) and predicted (asterisks = stochastic 
model, circles = deterministic model) spectra on September 10 at 6 instrument locations. 
The initial spectrum (H^ = 0.5 m) is indicated in each panel with a dashed line. 

Figure 8. Comparison of observed and predicted spectra on September 15 (H^ = 
0.4 m) (same format as Figure 7). 

Figure 9. Cross-shore evolution of the spectral levels at the peak frequency and 
the first four harmonics on September 15. The solid lines are the observed levels, the 
dotted lines are the deterniinistic model predictions, and the dashed lines are the stochastic 
model predictions. 

Figure 10. Bispectra observed (right panels) and predicted (left panels, from the 
stochastic model) offshore of the bar (sensor pl7, upper panels) and on the bar (sensor 
pl4, lower panels) on September 15 . The format of the panels is the same as in Figure 
2 . 

Figure 11. Bispectra observed (right panels) and predicted (left panels, from the 
stochastic model) on September 15 at sensors p4 (upper panels) and p3 (lower panels), 
both located on the beachface close to shore. The format of the panels is the same as in 
Figure 2. 

Figure 12. Comparison of observed and predicted spectra on September 24 (H^ 

= 0.8 m) (same format as Figure 7). 

Figure 13. Cross-shore evolution of the spectral levels at the peak frequency and 
the first two harmonics on September 24 (same format as Figure 9). 
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Figure 14. Bispectra observed (right panels) and predicted (left panels, from the 
stochastic model) on the bar (sensor 14, upper panels) and inshore of the bar (sensor p23, 
lower panels) on September 24 . The format of the panels is the same as in Figure 2. 

Figure 15. Comparison of observed and predicted spectra on September 21 
= 0.8 m) (same format as Figure 7). 
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